Different temporal trends in vascular plant and bryophyte communities along elevational gradients over four decades

Abstract Despite many studies showing biodiversity responses to warming, the generality of such responses across taxonomic groups remains unclear. Very few studies have tested for evidence of bryophyte community responses to warming, even though bryophytes are major contributors to diversity and functioning in many ecosystems. Here, we report an empirical study comparing long‐term change in bryophyte and vascular plant communities in two sites with contrasting long‐term warming trends, using “legacy” botanical records as a baseline for comparison with contemporary resurveys. We hypothesized that ecological changes would be greater in sites with a stronger warming trend and that vascular plant communities, with narrower climatic niches, would be more sensitive than bryophyte communities to climate warming. For each taxonomic group in each site, we quantified the magnitude of changes in species' distributions along the elevation gradient, species richness, and community composition. We found contrasted temporal changes in bryophyte vs. vascular plant communities, which only partially supported the warming hypothesis. In the area with a stronger warming trend, we found a significant increase in local diversity and dissimilarity (β‐diversity) for vascular plants, but not for bryophytes. Presence–absence data did not provide sufficient power to detect elevational shifts in species distributions. The patterns observed for bryophytes are in accordance with recent literature showing that local diversity can remain unchanged despite strong changes in composition. Regardless of whether one taxon is systematically more or less sensitive to environmental change than another, our results suggest that vascular plants cannot be used as a surrogate for bryophytes in terms of predicting the nature and magnitude of responses to warming. Thus, to assess overall biodiversity responses to global change, abundance data from different taxonomic groups and different community properties need to be synthesized.

Understanding variation among taxonomic groups in their responses to environmental change, including climate, is crucial for identifying priorities in conservation. However, the spatial variability of climate change complicates comparisons between studies of different taxa in different regions, and it is difficult to predict a priori differences between taxonomic groups. Thus, our knowledge is limited in identifying the taxonomic groups most sensitive to change in environmental conditions and, therefore, most in need of conservation action. We operationally define "sensitivity" here as responsiveness: that is, the degree to which a given community property changes in the face of change in environmental conditions.
For several reasons, we might expect vascular plants and bryophytes to show different responses to various sources of environmental change (Bagella, 2014;Becker-Scarpitta et al., 2017;Lee & La Roi, 1979;Möls et al., 2013;Vanneste et al., 2017). Compared with vascular plants, bryophytes are distinguished by their small size, high sensitivity to the moisture and chemistry of their immediate microenvironment (i.e., they are poikilohydric), lower temperature optima for growth, absence of roots and an efficient vascular system, and dispersal via spores (Glime, 2007). These two groups can show contrasted spatial patterns of diversity (Lalanne, Bardat, Lalanne-Amara, Gautrot, & Ponge, 2008;Mateo et al., 2016). For instance, vascular plants show a clear latitudinal diversity gradient of decreasing species richness with increasing latitude, while this is not true for bryophytes, for which temperate latitudes are equally diverse as tropical latitudes (Geffert, Frahm, Barthlott, & Mutke, 2013;Mateo et al., 2016). Some studies have also observed different patterns of community β-diversity (Kraft et al., 2011;Mateo et al., 2016) with vascular plant communities sometimes showing higher β-diversity along elevation gradients than bryophytes, suggesting a broader tolerance of bryophyte species to temperature (Glime, 2013;Lee & La Roi, 1979;Vanneste et al., 2017;Vittoz et al., 2010) or differences in dispersal ability (Lenoir et al., 2012).
General predictions for the effects of warming on vascular plants include declines in the abundance of cold-adapted species, an upward expansion of elevational range limits for warm-adapted species (Rumpf et al., 2018), and an increase of local species richness in area without water stress (Harrison, 2020;Steinbauer et al., 2018;Vellend, Dornelas, et al., 2017). Compared with vascular plants, some studies have suggested that changes in bryophyte communities are more strongly influenced by stochastic processes or by microenvironmental variation than macro-environmental conditions (Fenton & Bergeron, 2013;Pharo & Vitt, 2000;Raabe et al., 2010).
Because bryophytes have wider temperature tolerances and higher affinities to micro-environmental than to macro-environmental conditions, we might expect bryophytes to show lower sensitivity to warming than vascular plants. The consequences of warming for βdiversity are more difficult to predict given a paucity of studies on this topic (Socolar et al., 2016) (Nascimbene & Spitale, 2017). However, numerical simulations show that species with high dispersal capacity will be favored in areas experiencing strong environmental changes, in which case we might predict a decrease of β-diversity (Mouquet & Loreau, 2003) and thus biotic homogenization, at least temporarily (Clavel, Julliard, & Devictor, 2011).
As in many parts of the world, eastern Canada has shown a general warming trend over the past ~50 years, but with a strong eastwest gradient in the magnitude of warming in the province of Québec (Yagouti et al., 2008). Becker-Scarpitta, Vissault, and Vellend (2019) showed that the magnitude of temporal changes in vascular plant communities in three protected areas generally increased from east

T A X O N O M Y C L A S S I F I C A T I O N
Biodiversity ecology to west in southern Québec, with greater changes in community composition in areas of stronger warming in the last four decades.
For two of these three parks, both mountainous, the historical data also included bryophytes, thus presenting an opportunity to test for differential sensitivity among taxonomic groups to warming. Forillon National Park is located at the eastern tip of the province of Québec where warming has been less pronounced than in Mont-Mégantic Provincial Park, which is in central Québec where the warming trend has been greater (Yagouti et al., 2008; Supplementary material S1).
Parc Forillon has experienced an increase of 0.12°C ± 0.010°C/decade, while at Megantic the increase per decade is almost twice as strong: 0.20°C ± 0.014°C/decade.
Here, we report one of the first studies comparing long-term change in bryophytes and vascular plant communities in sites with contrasting long-term warming trends. In each of the two parks, we revisited ~50 legacy vegetation plots initially surveyed in the 1970s, applying the same methods as the original surveys. To minimize potentially confounding factors, plots were selected in mature forest ecosystems that have not experienced major natural or anthropogenic disturbances during the time between surveys. We test two main hypotheses: (i) For both bryophytes and vascular plants, the park with a stronger warming trend (Mont-Mégantic) has experienced greater long-term community changes than the park with a weaker warming trend (Forillon); (ii) vascular plant communities are more sensitive than bryophyte communities to climate warming. For each taxonomic group in each park, we quantified the magnitude of changes in (a) species' distributions along the elevation gradient, (b) species richness (alpha-diversity), and (c) community composition.
Note that these predictions were already tested for vascular plants in Becker-Scarpitta et al. (2019), to which we add data for bryophytes in this paper. and Acer saccharum Marsh. and Betula alleghaniensis Britt. at lower elevations (Majcen, 1981). Mont-Mégantic Provincial Park was created in 1994 (logging ceased in the 1960s during park planning) and covers ~55 km 2 . Our study plots span an elevational gradient from ~460 and 1100 m a.s.l., along which the vegetation gradient strongly resembles that in Parc Forillon: from temperate deciduous forests dominated by Acer saccharum., Fagus grandifolia Ehrh. and Betula alleghaniensis, to boreal forest dominated by Abies balsamea, and Picea rubens Sargent, Silva. (Marcotte & Grandtner, 1974). Given individual organisms that can live for decades, we expect considerable time lags between climate change and community-level responses. As such, we focus on climate change over a period of time offset from our vegetation surveys by 5-10 years , which aligns with a previous study showing a sharp warming gradient in the region (Yagouti et al., 2008). Average temperatures at Mont-Mégantic show an increase of ~0.90°C between 1960 and 2005 while at Forillon the estimated increase is about ~0.55°C for the same period (Supplementary Material S1).
In each plot, the authors listed all species in different strata (canopy trees, shrubs, herbs, and ground bryophytes) and each species was assigned an abundance coefficient following Braun-Blanquet et al. (1952). For vascular plants, our analyses focused on shrubs and herbs, which were combined into a single "understorey" stratum. All bryophyte species that were found on the ground were recorded (i.e., on organic litter and soil surface mineral layers, not including obvious deadwood, tree trunks and rocks); these surveys did not involve intensive searches for individual stems of rare species within moss carpets (i.e., some locally rare species were missed). After con- Original surveys in Forillon Park were conducted in June-September 1972 in 256 vegetation plots of 500 m 2 (Majcen, 1981).
We resurveyed a subset of 49 plots during July and August of 2015.
Original surveys at Mont-Mégantic were conducted in 1970 in 94 plots, roughly half of which fall outside the current park boundaries.
Plots were 400 m 2 in coniferous forest and 800 m 2 in broadleaved forests (Marcotte & Grandtner, 1974). We resurveyed a subset of

| Taxonomical database
Our taxonomical reference was the Taxonomic Name Resolution Service v4.0 (assessed in Feb 2017: http://tnrs.iplan tcoll abora tive.org) for vascular plants and Flore des bryophytes du Québec-Labrador (Faubert, 2012(Faubert, , 2013(Faubert, , 2014 for bryophytes. Our data set was collected by four different survey teams: one for each of the two original surveys: Forillon: (Majcen, 1981); Mont-Mégantic: (Marcotte & Grandtner, 1974); one for the recent Mont-Mégantic vascular plant survey: (Savage & Vellend, 2015); and one for the recent Mont-Mégantic bryophyte survey and for the recent survey of both taxonomic groups at Forillon (A. Becker-Scarpitta and assistants). Most plants were identified to the species level in the same way across surveys, so for these species the only harmonization step was to standardize names. Coarser levels of taxonomic resolution were used in some but not all surveys for certain species (e.g., a pair of similar species not identified to the species level), and for other species (e.g., spring ephemeral species) the timing of different surveys created doubt about the likelihood of comparable detection. In these situations, comparability was maximized by using the coarser level of resolution applied to all data sets, or by removing species (see Supplementary Material S2 for details on taxonomic standardization). All specimens identified at the species level were deposited in the Marie-Victorin herbarium (Université de Montréal, Canada.), and all locations were entered into the GBIF database (GBIF -https://www.gbif.org/).

| Statistical analysis
To minimize the effect of potential bias due to uncertain species identification or very rare species, all statistical analyses used presence-absence data for the subset of species with at least two occurrences in each time period for each park.
First, to test for upward elevational shifts in species distributions along the elevational gradients, we calculated the mean elevation across occurrences for each species (Z), in each park, at each time.
We modeled main effects and two-and three-way interactions between park (P), time (Y), and taxonomic group (T), each as a categorical variable. We thus tested for upward elevational shifts in species distributions, and how these shifts depend on park and taxonomic group. The model was given by: Where Z is the mean species elevation, β 1−8 are parameters to be estimated, b sp and b park are random intercepts accounting for the effect of species and park, and ϵ is the model residuals assumed to be normally distributed and independent. Following our first hypothesis, we expected the temporal change at Mégantic to be stronger than at Forillon (significant β 6 parameter). The second hypothesis predicted vascular plants to be more sensitive than bryophytes (significant β 7 parameter). All other terms of the model that are not part of our hypotheses ensure accurate estimation of other parameters in the model. Second, we used the exact same model structure and hypotheses to analyze temporal trends of αdiversity of both bryophytes and vascular plants in the two parks.
Third, we tested for temporal changes in community composition using two frameworks. In one framework, we first tested the multivariate homogeneity of group dispersions-that is, a multivariate β-diversity metric, a measure of dissimilarity-using the Jaccard dissimilarity index (PERMDISP, function betadisper, package "vegan" v.2.4-4; [Anderson, Ellingsen, & McArdle, 2006]). Then, we quantified the change in community composition over time using a permutational analysis of variance: PERMANOVA (Anderson, 2001) with Jaccard distances using 999 permutations to (functions adonis, package "vegan," and pairwise.adonis, package "pairwiseAdonis"). We first tested the significance of the three-way interaction (park * taxonomic group * time) in a full model (F = 4.45, p < .01), and second, we ran separate models for each combination of taxonomic group and park. The R 2 values computed from these separate PERMANOVA models were used to compare the magnitude of temporal change between groups. We used non-metric multidimensional scaling (nMDS) with Jaccard distances for visualization (function metaMDS, package "vegan").
In a second framework, to focus on the temporal component of change in species composition, we calculated for each plot the temporal pairwise β-diversity using the binary Jaccard dissimilarity index (for presence-absence data) between original and recent surveys (function beta.temp, package "betapart" v.1.5.2, Baselga & Orme, 2012). This framework allowed us to decompose the Jaccard dissimilarity into turnover (i.e., species replacement between plots) and nestedness (i.e., species loss between plots) (Baselga, 2010).
As a complementary analysis, we identified for each park separately the indicator species that characterized each period using the IndVal procedure (function multipatt, package "indicspecies"; De Cáceres et al., 2010). The function quantifies and tests the strength of the relationship between species occurrences and abundance and groups of sites. To determine whether a species is associated with a group of sites, a statistical test should be performed to evaluate the null hypothesis that this relationship does not exist. The null hypothesis, a species that was observed at a site belonging to the target site group, is only due to chance. All statistical analyses were performed in R v.3.4.2 (R Foundation for Statistical Computing 2018).

| Species distributions along elevation gradients
Among the four taxonomic group-site combinations, we failed to detect any significant upward shifts in mean species elevations over time-that is, none of the main effects or interaction terms including time were significant (Table 1a)  Despite the lack of trends in mean species elevation for bryophytes, there was substantial variation among species, more than for vascular plants (Table 1a,

| Temporal changes in species richness
At Forillon, the total number of bryophyte species across plots was greater in the recent survey (38 species) than in the original survey We detected a significant increase in mean plot-level species richness (α-diversity) only for vascular plants at Mont-Mégantic (Table 1b). For bryophytes in both parks, and for vascular plants at Forillon α-diversity showed no changes over time (Table 1b).

| Temporal shift in community composition and heterogeneity
We found significant temporal changes in vascular plant β-diversity and community composition at Mont-Mégantic (   Table 1a for statistical analysis. Note: β-diversity was measured as the multivariate distance of plots to the time-specific centroid in multivariate space. The R 2 of the shift in community composition reflects the proportion of variation in community composition explained by time. Jaccard temporal β-diversity values are decomposed into turnover and nestedness between original and recent surveys. Bold values indicate significant statistical differences p < .05.

TA B L E 2 Temporal changes in β-diversity (dissimilarity) and community composition
For both sites, we found significant temporal shifts in bryophyte community composition; for vascular plants, there was a significant compositional change only at Mont-Mégantic (Table 2, Figure 2). As predicted, the magnitude of the vascular plant compositional shift was greater for Mont-Mégantic (R 2 = 0.037) than Forillon (R 2 = 0.019).
However, bryophyte communities experienced an equal magnitude of compositional shift for both Forillon and Mont-Mégantic (R 2 = 0.07 and 0.074, respectively). Moreover, compositional shifts were greater for bryophytes than for vascular plants at both sites ( Table 2).
Temporal β-diversity showed contrasted patterns between taxonomic groups. Jaccard dissimilarity between original and recent surveys was higher for bryophytes than for vascular plants. The greater change in temporal β-diversity for bryophytes in the two parks is consistent with the PERMANOVA results for changes in community composition. For all combinations of taxonomic groups and parks, species turnover contributed more than nestedness to overall changes. However, nestedness was considerably higher for vascular plants than for bryophytes ( Table 2). In support of the β-diversity analysis, Table 3 shows the indicator species that characterized each

F I G U R E 2 Temporal changes in (a) bryophyte and (b) vascular plants community composition over time in Forillon and Mont-
Mégantic, based on Jaccard's dissimilarity and non-metric multidimensional scaling (NMDS). At both sites, bryophytes showed significant shifts in composition over time and a significant increase in β-diversity. For vascular plants, significant shifts in community composition and a significant decrease in β-diversity were found at Mont Mégantic but not Forillon ( Table 2).   (Table 1b), and a stronger change in β-diversity and community composition for vascular plans as well (Table 2).
However, we did not detect any shift in species distributions along the elevation gradient (Table 1a, Figure 1) and the higher magnitude of changes in bryophyte community composition at both sites was not expected (Table 2, Figure 2).

| Patterns along the gradient of warming trends
At Mont-Mégantic, where the warming trend has been strongest, we found a significant increase in α-diversity (Table 1b), a significant decrease in β-diversity, and a stronger shift in community composition than in Forillon (Table 2, Figure 2). At Forillon, where the regional warming trend has been weaker, we found neither shifts in elevational distributions (Table 1a, Figure 1) nor changes in αdiversity for either vascular plants or bryophytes (Table 1a). The lack of upward shift in elevation of vascular plants in response to warming is clearly explained by the use of presence/absence data and contrasts with our previous results using abundance data (Becker-Scarpitta et al., 2019) in addition to studies from other sites (Bertrand et al., 2011;Chen et al., 2011;Lenoir et al., 2008). The temporal increase in diversity of vascular plants at Mont-Mégantic is also consistent with the prediction that warming should lead to increased local diversity in areas without severe moisture stress (Harrison, 2020;. Our results suggest that a temperature increase of ~0.9°C (over the last 45 years at Mont-Mégantic) does not have as strong an impact on the local diversity of bryophytes as it does for vascular plants, assuming that the richness increase at Mégantic was indeed due to warming. This interpretation is also supported by the absence of any relationship between bryophyte richness and elevation in the two parks (Forillon estimate = 0.0042 ± (SE) 0.0024, t = 1.8, p = .079; Mégantic estimate = 0.0022 ± (SE) 0.0014, t = 1.6, p = .11) and consistent with previous results in the literature (Bruun et al., 2006, Grytnes, Heegaard, & Ihlen, 2006, Odland et al., 2014but see Vittoz et al., 2010). The decrease in β-diversity over time for vascular plants at Mont-Mégantic is consistent with the hypothesis that warming might cause biotic homogenization (Socolar et al., 2016;Urban, 2015).
Note that this result using presence-absence data was different to the result in Becker-Scarpitta et al. (2019) using untransformed abundance data (for which we found no homogenization), but similar to the finding of homogenization using fourth-root transformed abundances in Savage and Vellend (2015). These results collectively illustrate that the locally dominant species-whose influence is minimized or eliminated via fourth-root or presence-absence transformationcan mask homogenization created by species of lower abundance.
This result highlights how ecological conclusions are dependent on the data transformation, even with widely used diversity metrics.

| Sensitivity of bryophytes vs. vascular plants
We

| Potential non-climatic drivers of vegetation change
Our sites were chosen specifically due to their contrasting warming trends and lack of other obvious major drivers of vegetation change, but there are certainly other possible drivers of ecological change that might play a role in this region. Among the indicator species of recent surveys (Table 3), two were non-native: Galeopsis tetrahit and Epipactis helleborine, the latter of which has increased considerably in recent decades throughout its North American range, even in western Canada (Marie-Victorin, 1997;. As such, some vegetation changes might be due simply to protracted periods of non-native species expansions, regardless of local environmental change. Another potential factor is changes in white-tailed deer browsing, which has increased over the past century in much of North America (Côté, Rooney, Tremblay, Dussault, & Waller, 2004). The indicator species of the recent survey at Mont-Mégantic (Table 3) include species known to benefit from high levels of deer browsing: Dennstaedtia punctilobula and Carex spp. (Augustine & Decalesta, 2003;de la Cretaz & Kelty, 2002;Frerker, Sabo, & Waller, 2014;Rooney, 2009). Interestingly, at Forillon, deer are actually thought to have decreased in abundance due to the expansion of the coyote population in the 1970s (UQCN, 2005), and we found no such species associated with recent surveys in Forillon Park.
Our most difficult result to interpret was the strong species turnover of bryophyte communities at Forillon, which has not experienced strong long-term increase in temperature, precipitation, or atmospheric nutrient deposition (Commission Joint International, 2014;Hember, 2018). We can only speculate regarding potential non-climatic drivers of bryophyte community change.
First, as in all legacy studies, there is the potential for observer biases due to (i) different sampling effort between original and recent surveys, or (ii) species' identification errors. It seems highly likely that detection probabilities and the potential for identification errors are greater for bryophytes than for vascular plants, although we have no reason to suspect this caused systematic increases or decreases in particular species frequencies (necessary to explain overall compositional shifts). We paid very close attention to repeating the original survey methods, focusing on the visually obvious species in each microsite (i.e., not examining each individual moss stem closely on the field), and the lack of any difference over time in local richness suggests comparable species' detection abilities in the two surveys. Although we cannot exclude the possibility that a real richness change was canceled out by a change in survey effort, we have no reason to suspect this rather unlikely coincidence. Second, given uncertainty in the comparability of abundance estimates across time for bryophytes, we also decided to use presence-absence data. In short, changes in observer effort seem highly unlikely to account for the temporal change in species composition. Furthermore, our procedure to homogenize taxonomy across data sets was quite conservative to reduce bias due to misidentifications.
One potential hypothesis to explain compositional change over time in bryophyte communities relates to the history of park protection. Forillon Park was established (and so protected) only 2 years before the original survey was conducted, and parts of the park previously included homesteads (i.e., non-forest land uses).
Mont-Mégantic was established as a park more recently (1994), but logging activities (the only prominent land use) ceased ~15 years before the original survey. Although plot selection focused only on non-disturbed forests, delayed responses to historical disturbance or metacommunity dynamics involving dispersal of species from sites undergoing rapid succession may have caused local shifts in composition and increased in dissimilarity (β-diversity).
It has been previously shown that managed forests tend to have a lower β-diversity than protected forests (Kaufmann, Hauck, & Leuschner, 2017). The increase in bryophyte β-diversity might partially be due to the recovery of natural forest that occurring in the 1970s.
There is also the possibility that changes in bryophyte communities resulted from interactions with changing vascular plant communities. Studies have shown that bryophyte diversity and abundance can be negatively correlated with total vascular plant biomass (Virtanen, Eskelinen, & Harrison, 2017), cover (Jiang, Liu, Song, Yu, & Shao, 2015)

| Conclusion and conservation implications
Overall, we found a significant temporal shift in the composition of both taxonomic groups in both parks but only one significant change in α-diversity. This study shows the limitation of using presenceabsence data to detect shifts in distribution when using a diachronic approach with only two time points. While abundance data provided the necessary power to detect changes in species distributions along this elevational gradient (Becker-Scarpitta et al., 2019), presence/absence data clearly did not. Our results are in accordance with recent meta-analyses and syntheses showing that local diversity can remain unchanged (or increase or decrease with equal likelihood) despite strong changes in composition (Blowes et al., 2019;Dornelas et al., 2014;Magurran et al., 2018;Spaak et al., 2017;Vellend, Dornelas, et al., 2017). Finally, regardless of whether one taxonomic group is systematically more or less sensitive to climate change than another, our results suggest that one taxonomic group (e.g., vascular plants) cannot be used as a surrogate for others (e.g., bryophytes) in terms of predicting the nature and magnitude of responses to climate change (Bagella, 2014;Becker-Scarpitta et al., 2017;Outhwaite et al., 2020). In the same plots that experienced the same climate change, we found that communities of bryophytes and vascular plants did not predictably change in the same ways (Bagella, 2014;Becker-Scarpitta et al., 2017;Lalanne et al., 2008;Odland et al., 2014;Slack, 1977).
Thus, to assess overall biodiversity responses to global change, abundance data from different taxonomical groups and different community properties need to be synthesized.

ACK N OWLED G M ENTS
We would like to thank the field assistants who greatly contributed to the collection of data in the field, the identification of species in the laboratory, and the realization of the herbarium: Melissa Paquette and Sara Gaignard. This work was made possible thanks to the support of park managers, especially Camille-Antoine Ouimet (at Mont-Mégantic) and Daniel Sigouin (at Forillon). We also thank Guillaume Blanchet, Michael Belluau, and Willian Vieira for constructive discussion on this project. The project was funded by the Natural Sciences and Engineering Research Council, Canada, and field mission was supported by the Quebec Center for Biodiversity Sciences (QCBS). writing -original draft (equal); writing -review and editing (equal).

CO N FLI C T O F I NTE R E S T
The authors state that they have no conflicting interests.